rm(list = ls())
here::i_am(file.path("code", "13_fs_alternatives.R"))
library(here)
source(here("code", "config.R"))

est_df <-  read_parquet(here("data", "analysis", "analysis.parquet"))

ols <- feols(
  is_naacp ~ vet_combined |
    birthyr_cards + bpl_cards + board_identifier +
    exemption^married_cards + farmer + laborer + farm_laborer,
  data = est_df,
  vcov = "hetero"
)

rhs <- model.matrix(ols)
lhs <- model.matrix(ols, type = "lhs")
my_nonvet <- mean(lhs[rhs == 0])

ols_coef <- sprintf(
  "%.04f",
  coefficients(ols)
)

ols_t <- sprintf(
  "%.02f",
  tstat(ols)
)

res1 <- feols(
  is_naacp ~ 1 | board_identifier + birthyr_cards + bpl_cards + exemption^married_cards +
    farmer + laborer + farm_laborer |
    vet_combined ~ order_num_serial_norm,
  data = est_df,
  cluster = ~ serial_num
)

res2 <- feols(
  is_naacp ~ 1 | board_identifier + birthyr_cards + bpl_cards + exemption^married_cards +
    farmer + laborer + farm_laborer |
    vet_combined ~ order_num_norm,
  data = est_df |>
    mutate(
      order_num_norm = local_order_num/registrants_second_report,
    ),
  cluster = ~ serial_num
)

res3 <- feols(
  is_naacp ~ 1 | board_identifier + birthyr_cards + bpl_cards + exemption^married_cards +
    farmer + laborer + farm_laborer |
    vet_combined ~ local_order_num_from_serial,
  data = est_df,
  cluster = ~ serial_num
)

res4 <- feols(
  is_naacp ~ 1 | board_identifier + birthyr_cards + bpl_cards + exemption^married_cards +
    farmer + laborer + farm_laborer |
    vet_combined ~ local_order_serial_rank_norm,
  data = est_df,
  cluster = ~ serial_num
)

etable(
  res1,
  res2,
  res3,
  res4,
  tex = TRUE,
  signif.code = c("***" = 0.01, "**" = 0.05, "*" = 0.10),
  fitstat = ~ n + r2 + kpr,
  replace = TRUE,
  file = file.path(tab_dir, "iv_robustness_instruments.tex"),
  fixef.group =  list(
    "-_Birth year and state" = c("Birth year", "Birth state"),
    "-_Exemption claim $\\times$ married" = c("Exemption-Married"),
    "-_Prewar occupation" = c("Farmer", "Laborer", "Farm laborer")
  ),
  headers = c(
    "Baseline",
    "Order num.",
    "Not scaled",
    "Rank"
  ),
  extralines = list(
    "__Dep. var. mean (nonveterans)" = rep(my_nonvet, 4),
    "__OLS coefficient" = rep(ols_coef, 4),
    "__OLS $t$-statistic" = rep(ols_t, 4)
  ),
  postprocess.tex = function(x) { 
    x <- gsub("Kleibergen-Paap", "First stage $F$-statistic", x, fixed = TRUE)
    return(x)
  }
)
